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This  paper  presents  the  results  of  flow  computations  around  plunging  airfoils  using  a 
low  dissipation  finite  volume  scheme  that  preserves  kinetic  energy.  The  kinetic  energy  pre¬ 
serving  scheme  is  briefly  described  and  numerical  experiments  are  done  on  low  Mach,  low 
Reynolds  numbers  plunging  airfoils.  Results  are  discussed  and  compared  with  experimental 
data. 


Nomenclature 


Re  Reynolds  Number 

Pr  Prandtl  Number 

M  Mach  Number 

Sr  Strouhal  Number,  ujhL/V, ^ 
x  Spatial  Coordinate 

t  Time  variable 

v*  1  ith  component  of  flow  velocity 

Voo  Freestream  velocity 

p  Density 

p  Pressure 

E  Total  Energy 

H  Specific  Total  Enthalpy 

k  Specific  Kinetic  Energy 

p  Dynamic  Viscosity 

A  Volume  Viscosity 


k  Coefficient  of  heat  conduction 

Cp  Specific  Heat 

L  Chord  Length 

h  Nondimensional  Plunge  Amplitude 

(j  Angular  frequency 

k  Reduced  frequency,  loL/Voq 

Q  Fixed  spatial  domain 

dQ  Boundary  of  the  domain  Q 

Superscript 

i,j  Coordinate  directions 

Subscript 

o,p  Grid  Locations 

oc  Far  field  quantity 


I.  Introduction 

Over  the  last  century,  steady  airfoils  have  been  thoroughly  studied  yielding  important  experimental  and 
numerical  results.  On  the  contrary,  unsteady  flows  around  airfoils  have  not  been  well  characterized  due  to 
their  higher  complexity.  It  is  still  a  challenge  to  obtain  accurate  empirical  data  and  current  computational 
capabilities  are  not  sufficient  to  gain  high  resolution  of  the  phenomenon.  In  recent  years,  unsteady  airfoils 
have  been  gaining  a  lot  of  interest,  especially  towards  the  use  of  flapping  flight  in  the  development  of  micro 
unmanned  aircraft  vehicles  (UAVs). 

The  flow  around  a  plunging  airfoil  is  characterized  by  the  generation  of  vortices  that  strongly  interact 
in  the  wake  of  the  airfoil.  In  order  to  study  such  an  unsteady  flow  numerically  and  capture  the  phenomena 
present  in  the  wake  with  good  resolution,  it  is  crucial  to  use  a  scheme  that  introduces  as  little  dissipation 
as  possible.  When  the  artificial  dissipation  of  the  scheme  becomes  too  large,  there  is  significant  uncertainty 
whether  the  damping  observed  is  due  to  natural  viscosity  or  numerical  dissipation. 

Jameson  recently  proposed  a  scheme  that  could  overcome  this  problem.  In  his  Kinetic  Energy  Preserving 
(KEP)  scheme,  the  global  balance  of  kinetic  energy  over  the  computational  domain,  along  with  the  balance 
of  mass,  momentum  and  total  energy  is  respected.  The  KEP  scheme  provides  improved  stability,  removing 
almost  completely  the  need  of  adding  artificial  viscosity  and  making  it  the  best  candidate  for  this  study. 

In  this  work,  we  used  the  KEP  scheme  to  perform  computations  of  flows  around  plunging  airfoils  at  low 
Mach  and  low  Reynolds  number.  Section  2  of  this  paper  briefly  describes  the  Kinetic  Energy  Preserving 
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scheme  developed  by  Jameson  and  the  formalism  used  to  solve  the  equations.  Section  3  describes  how  the 
computations  were  performed.  In  Section  4,  results  are  presented  and  discussed. 


II.  Kinetic  Energy  Preserving  scheme  for  viscous  flow 

In  a  recent  paper,  Jameson1,2  presented  a  new  semi-discrete  finite  volume  scheme  to  solve  the  Navier- 
Stokes  equations.  This  scheme  has  the  property  to  satisfy  the  global  conservation  law  for  kinetic  energy.  We 
shall  briefly  describe  this  scheme  in  the  present  section. 


II. A.  Continuous  Model 


First,  consider  the  three-dimensional  Navier- Stokes  equations  in  their  conservative  form: 


where 


du 

dt 


_9_ 

dxi 


f(u)  =  0 


’  p  " 

pv1 

pv1 

pvlv1  +  pS%1  —  a11 

pv2 

and  fl  — 

pvlv 2  +  pSl2  —  Gl2 

pv 

pvlv 3  +  p5z3  —  a1 3 

-PE. 

pvzH  —  vi  <jli  —  qi  _ 

(1) 


(2) 


The  viscous  stress  tensor  is  given  for  a  Newtonian  fluid  by  =  A +  p  (fyj  +  § Often 

in  aerodynamics,  A  is  taken  to  be  equal  to  —  |/i.  The  heat  flux  is  proportional  to  the  temperature  gradient 
(Fourier’s  Law)  qi  =  — 

An  equation  for  the  kinetic  energy  k  =  ^pvl<2  can  be  derived  by  combining  the  continuity  and  the 
momentum  equations.  Indeed, 


dk 

dt 


d_ 

dt 


( pyi )  - 


vl2  dp 

~Y~dt 


It  follows  by  substituting  (pv2)  and  ^  by  their  corresponding  fluxes  that: 


dk  d 
dt  +  dxi 


dvi  •  •  dvl 

V - r  —  a  J - 7 

dxi  dxi 


(3) 


We  assume  that  we  are  interested  in  a  domain  Q  fixed  in  space,  dfl  denotes  the  boundary  of  fh  By 
integrating  (3)  over  the  domain  we  get  a  global  conservation  law  for  kinetic  energy. 


Q  dQ 


v  a 


njdS  + 


dV 


(4) 


Definition  1  A  numerical  scheme  to  solve  the  viscous  Navier-Stokes  equations  is  said  to  be  Kinetic  Energy 
Preserving  if  it  satisfies  a  discrete  analog  of  (4 )■ 

Here  we  have  assumed  that  the  domain  contains  no  discontinuity.  If  a  shockwave  is  present  in  SJ,  the 
relation  (4)  does  not  hold  anymore. 
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II. B.  Semi-discrete  approach 


Now,  we  consider  a  finite  volume  discretization  of  the  governing  equations  in  the  domain  It  .  The  generic 
cell  is  a  polyhedral  control  volume  o.  Each  cell  has  one  or  more  neighbors.  The  face  separating  cell  o  and 
cell  p  has  an  area  Aop ,  and  we  define  nlop  to  be  the  unit  normal  to  this  face,  directed  from  o  to  p.  Evidently 
nlop  =  —npo.  We  also  define  Slop  =  Aopnlop.  S%op  can  be  interpreted  as  the  projected  face  area  in  the  coordinate 
direction  i. 

Boundary  control  volumes  are  closed  by  an  outer  face  of  directed  area  Sl0  —  —  Slop  (a  control  volume 
is  delimited  by  a  closed  surface). 

In  this  framework,  the  semi-discrete  finite  volume  approximation  of  the  governing  equations  takes  the 
form: 

vol°ijt+  E  fop-niPAop  =  o  (5) 

p  neighbor 

or 

vol°^r+  E  fop  ■  Sip  =  0  (6) 

p  neighbor 

For  a  boundary  control  volume  6,  another  contribution  to  the  fluxes  /£  •  S&  comes  from  the  outer  face. 

Now  we  assume  that  uQ  and  flop  take  the  form: 


po 

1 

a, 

o 

E 

PoVl 

(pv'v^op  +  (p5a  -  <ja)op 

PoV20 

and  fop  = 

(/ pvlv2)op  +  (p5l2  -  al2)op 

PoVl 

(, pvlv3)op  +  (; pS* 3  -  at3)op 

-PoEo. 

.  (pvlH)op  -  (v3 g13  +  q3)op _ 

Jameson  has  exhibited  a  set  of  sufficient  conditions  on  the  elements  of  f%  that  lead  to  a  Kinetic  Energy 
Preserving  (KEP)  scheme. 

Proposition  1  If  the  elements  of  flop  defined  in  (15)  satisfy  the  following  conditions: 
a  -  (/wV)op  =  +  vi) 

b  -  (p513  -  crl3)op  =  |( pS 13  -  o%3)0  +  \{pb13  -  crt3)p 

and  if  the  fluxes  at  the  boundaries  are  evaluated  such  that: 
c  -  fl  =  f(uh)  where  b  is  a  boundary  control  volume 

then  the  semi  discrete  finite  volume  scheme  (6)  satifies  the  discrete  global  variation  law  for  kinetic  energy. 


Indeed  in  that  case,  the  discrete  kinetic  energy  kQ  satisfies  the  following  relation: 


Jt  E  VOl°k°  =  ~  E  Sb  (  vb  (  Pb  +  Pb\  )  -  <°b 


ST  (  v ^vI  +  vp  pi  i.7  vo  +  vp  pi 

/  >  I  P°  o  ^op  ao  2-^  9 


(8) 


which  is  indeed  a  discretization  of  (4). 
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Condition  a  of  the  previous  proposition  is  not  very  restrictive  and  allows  some  degrees  of  freedom  in  the 
construction  of  the  fluxes  defined  in  (15). 

Let’s  denote  by  gop  the  arithmetic  average  of  the  quantity  g  between  cell  o  and  cell  p  :  gop  =  (gQ  +  gp)/2. 
We  can  rewrite  condition  a 

(Pvivj)0p=  (pvi)opV3oP  (9) 

We  can  evaluate  the  average  (pvl)  by  any  means  (  (pvl)  =  Poptfop  or  Pyi  op  f°r  example)  and  then 
deduce  by  using  (9)  to  satisfy  condition  a. 

This  degree  of  freedom  could  be  used  to  design  schemes  that  also  satisfy  other  properties  than  Kinetic 
Energy  Preservation. 


III.  Direct  Numerical  Simulation  of  a  plunging  airfoil 

The  Kinetic  Energy  Preserving  scheme  described  above  was  used  to  compute  the  flow  around  a  plunging 
airfoil.  Computations  were  done  for  a  NACA  0012  airfoil  oscillating  in  a  uniform  flow.  The  transversal 
motion  of  the  airfoil  is  given  by  h(t)  =  h  •  cos(ut). 

The  freestream  flow  is  characterized  by  a  Mach  number  M^,  a  far  field  temperature  T ^  and  a  far  field 
density  p^. 

Viscosity  is  evaluated  using  Sutherland’s  formula 

7^3/2 

»m=c—s 

For  air,  at  reasonable  temperatures,  C  =  1.456  x  10-6kg/(  ms\/K)  and  S  =  110. 4K. 

The  Reynolds  number  is  based  on  the  chord  length  of  the  airfoil  L  and  the  free  stream  velocity  =  ^oocoo 

Re  =  f)o°LVo°  where  =  giT^)  (10) 


The  Prandtl  number  is  given  by 

Pr=^v  (11) 

Kj 

it  was  taken  to  be  equal  to  0.75. 


Computational  Domain 

Simulations  were  done  on  a  structured  UC-  mesh”  counting  4096  x  512  cells.  The  computational  domain 
extends  roughly  30  chord  lengths  downstream  and  20  chord  lengths  upstream,  as  can  be  see  on  figure  2. 
The  mesh  is  subject  to  rigid  body  motion  and  moves  with  the  airfoil. 


Numerical  fluxes 


The  convective  fluxes  have  to  be  modified  to  account  for  the  motion  of  the  mesh.  First,  let  us  consider  the 
hyperbolic  system  of  equation  (1)  =  0  and  integrate  it  over  the  moving  domain  Q  (£).  We  have 

using  the  divergence  theorem 

/  iftdV+  J  =  °  (12) 

Q(t) 


However,  since 


/ 


l/“iv+  / 

Q(t)  dQ(t) 


uvl  •  nldS 


(13) 
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(a)  Global  domain  (b)  Mesh  detail  near  the  trailing  edge 

Figure  1.  Computational  domain 


where  vl  is  the  speed  of  the  boundary  of  the  domain,  it  follows  that 

J  udV  +  J  ( fl(u )  —  uvz)  •  nldS  =  0  (14) 

Q(t)  dQ(t) 

Denote  vlop  the  velocity  of  the  edge  separating  cells  o  and  p.  The  convective  fluxes,  are  now  given  by 


f 

J  c 


op  convective 


'-’op  )  ]  op 

I'll  _1_  rr,  VI 


[p(v 

[p(v*  -<p)v1)\op+Pop& 
[p(vi  -Kp)v2)]op+popSi2 

[p(vi  ~Kp)v3)]op+PopSi3 

[p{vl  -vip)E+pv‘ 


J  op 


We  then  used  the  following  averaging  formula: 


^op)\op  Pop  (V op  Vop) 

[p(vl  -  «op)vJ]op  =  Pop(vlop  -  Kp)vJop 

\p(v  Kp  )  ^  ^  ]  Qp  Pop(vop  Kp  )  ^  op  Pop^'op 


(15) 


(16) 


Condition  a  of  proposition  (1)  does  not  require  a  specific  form  for  the  energy  flux.  We  defined  it  in  a 
consistent  manner  with  the  continuity  and  momentum  fluxes. 

Viscous  stress  was  evaluated  in  each  cell  by  introducing  a  complementary  mesh,  for  which  cell  vertices  are 
the  centers  of  the  original  control  volumes. 


Time  integration 


Time  integration  was  done  using  a  TVD  Runge  Kutta  second  order  multistage  time  stepping  scheme3.  For 
a  semi  discrete  law  in  the  form 


—  +  R(u,t)  =0 


(17) 
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the  scheme  advances  from  time  n  to  time  n  +  1  by 

u1  =un  -  A tR(un,tn) 
un+ 1  =  X- un  +  1 u 1  - 

The  explicit  dependance  in  time  of  the  operator  R  is  due  to  mesh  motion.  This  means  that  both  the 
location  and  velocity  of  the  mesh  need  to  be  updated  at  £n+1  before  evaluating  R{v}  Rn+1). 

This  time  discretization  does  not  guarantee  the  preservation  of  kinetic  energy  in  time.  One  could  use  a 
Crank-Nicholson  semi  implicit  scheme  as  suggested  by  Jameson1  to  ensure  conservation  in  time,  but  the 
computational  costs  would  increase  enormously. 


Far  field  artificial  dissipation 

As  the  mesh  coarsens  in  the  far  field,  roughly  5  chord  lengths  away  from  the  airfoil,  small  spurious  oscillations 
associated  with  acoustic  waves  can  be  observed.  It  was  shown  in  previous  work1,2,4  that  the  number  of 
cells  required  to  ensure  a  non  oscillatory  solution  and  stability  was  governed  by  the  local  cell  Reynolds 
number,  which  has  to  be  of  the  order  of  unity  to  guarantee  these  properties.  However,  covering  the  entire 
computational  domain  with  cells  as  fine  as  the  ones  in  the  wake  is  currently  too  expensive. 

A  small  amount  of  dissipation,  based  on  the  Jameson- Schmidt- Turkel  (JST)  scheme7,8  was  added  in 
the  far  field  to  control  these  unphysical  oscillations  and  prevent  the  explosion  of  the  number  of  grid  points. 
Furthermore,  at  the  very  edges  of  the  computational  domain,  where  the  cells  are  the  biggest,  artificial 
dissipation  is  larger  and  behaves  like  a  sponge  that  prevents  the  reflection  of  the  acoustic  waves  into  the 
computational  domain. 

The  dissipation  introduced  in  the  far  field  was  derived  from  the  JST  scheme  by  dropping  the  lower  order 
diffusive  term  and  conserving  the  higher  order  term  to  control  odd/even  modes.  It  shall  be  noted  that  no 
dissipation  was  introduced  in  the  near  field  of  the  airfoil,  an  area  encompassing  about  70%  of  the  cells.  If  we 
consider  the  conservation  equation  ^  =  0,  the  truncation  error  introduced  by  our  second  order 

scheme  can  be  seen  as  a  continuous  term  in  a  modified  differential  equation 

§2  +  ^/M  =  0(Ax»,Ata) 

The  idea  is  to  introduce  an  extra  diffusive  term  that  will  modify  the  truncation  error 


with  >  0  and  p  >  2  to  preserve  the  order  of  the  scheme.  If 


du  1 
dt  Ax 


K 


=  0 


is  a  finite  volume  semi-discretization  of  the  equation  where  hi±  i  is  the  numerical  flux,  we  can  introduce  a 
correction  di±  i  to  the  flux  to  obtain  the  desired  property.  This  can  be  done  by  taking 


di+i  =  ai+ie^  (~ui+ 2  +  3ui+1  —  3 iq  +  i) 

similarly  to  what  is  done  in  the  JST  scheme.  ai+ 1  is  proportional  to  the  spectral  radius  of  the  local  jacobian 
matrix  and  e^4)  is  a  switch  to  add  dissipation  only  where  needed,  as  described  in  the  original  JST  scheme. 
di+ 1  is  proportional  to  Ax3|^||.. 


IV.  Results 

We  now  present  the  results  obtained  using  the  numerical  methods  described  above.  We  chose  the  various 
parameters  describing  the  motion  of  the  plunging  airfoil  to  match  the  ones  studied  by  Jones  and  Platzer5,6. 
Flows  around  plunging  airfoils  can  be  classified  according  to  their  Strouhal  numbers  Sr  =  =  kk.  In  all 

the  results  presented,  the  Mach  number  is  taken  to  be  M ^  =  .2  and  the  Reynolds  number  is  Re  =  1850. 
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30 
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X-axis 


Figure  2.  Artificial  dissipation  is  added  only  in  the  darker  area 


Drag  production  at  low  Strouhal  number  -  Sr  =  0.29 

The  amplitude  of  the  plunging  motion  is  h  =  0.08  and  the  reduced  frequency  k  =  3.6  resulting  in  a  Strouhal 
number  Sr  =  0.288.  For  such  a  Strouhal  number  the  resulting  flow  exerts  drag  on  the  airfoil.  Figure  3 
shows  the  density  contour:  not  only  the  flow  pattern  is  clear  on  this  picture,  but  this  is  also  evidence  that 
there  are  some  non  negligeable  compressibility  effects,  even  at  such  a  low  Mach  number.  Figure  4  represents 
the  vorticity  of  the  flow.  As  can  be  seen,  the  vortical  structure  of  the  wake  is  in  agreement  with  the  results 
obtained  by  Jones  et  ah,  cf.  figure  5.  Figure  6  is  a  plot  of  the  evolution  in  time  of  the  lift  coefficient  Ci  and 
the  drag  coefficient  Cd-  A  positive  Cd  indicates  that  the  fluid  exerts  drag  on  the  airfoil  in  the  flow  direction. 
On  the  abscissa,  the  non  dimensional  time  is  defined  by  t  =  ^ ,  to  =  y~  . 


X  axis 


Density 


Figure  3.  Density  field  in  the  airfoil’s  wake  -  Sr  =  0.29 


Thrust  generation  at  high  Strouhal  number  -  Sr  =  0.60 

As  the  Strouhal  number  is  increased  the  flow  pattern  changes  and  thrust  is  generated.  In  this  section,  we 
present  the  results  obtained  for  a  Strouhal  number  of  Sr  =  0.60.  Figures  7  and  8  depict  the  computed  density 
and  vorticity  fields  for  h  =  0.1  and  k  =  0.6.  The  flow  pattern  is  in  excellent  agreement  with  experimental 
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Figure  4.  Vorticity  distribution  in  the  airfoil’s  wake  -  Sr  =  0.29 


Figure  5.  Streak  lines,  experimental  data  by  Jones  and  Platzer  -  Sr  =  0.29 


Figure  6.  Lift  and  Drag  history  -  Sr  —  0.29 


"O 


O 


data  by  Jones  et  al.  corresponding  to  this  Strouhal  number  (see  figure  9)  and  we  observe  on  figure  10  that 
thrust  is  indeed  generated.  However,  it  should  be  noted  that  Jones’  experimental  result  was  obtained  for 
h  =  0.2  and  k  =  0.3.  When  we  tried  to  compute  the  flow  for  these  values,  we  still  observed  generation  of 
thrust,  but  the  flow  pattern  was  quite  different.  The  large  amplitude  of  the  motion  (for  h  =  0.2)  creates 
important  leading  edge  vortices  that  interact  strongly  with  the  trailing  edge  vortices  in  the  wake,  as  can  be 
seen  on  figure  11. 
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Figure  7.  Density  field  in  the  airfoil’s  wake  -  Sr  =  .60,  h  —  .  1,  k  =  .6 
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Figure  8.  Vorticity  distribution  in  the  airfoil’s  wake  -  Sr  =  .60,  h  =  .1,  k  —  .6 


Figure  9.  Streak  lines,  experimental  data  by  Jones  and  Platzer  -  Sr  =  0.60,  h  —  .2,  k  =  .3 


Lift  and  thrust  generation  -  Sr  =  1.5 

This  is  by  far  the  most  interesting  case,  as  lift  and  thrust  are  generated  by  the  oscillating  airfoil.  The 
nondimensional  plunging  amplitude  is  h  =  .12  and  the  reduced  frequency  k  =  12.3,  resulting  in  a  Strouhal 
number  Sr  =  1.48.  Figures  12  and  13  depict  the  density  and  the  vorticity  of  the  flow  in  the  airfoil’s  wake. 
The  dual-mode  vortex  street  described  by  Jones  et  al.5  is  clearly  visible.  Our  numerical  computations  are 
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Time  history  of  Cd  and  Cl 


Figure  10.  Lift  and  Drag  history  -  Sr  =  0.60,  h  =  .  1,  k  —  .6 
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Figure  11.  Vorticity  field  in  the  airfoil’s  wake  -  Sr  —  .60,  h  =  .2,  k  =  .3 


again  in  excellent  agreement  with  Jones  et  al.  experimental  results,  as  shown  in  figure  14. 
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Figure  12.  Density  field  in  the  airfoil’s  wake  -  Sr  =  1.5 


Density 


V.  Conclusion 

This  paper  shows  how  to  solve  with  high  resolution  the  flow  around  plunging  airfoils  using  a  finite  volume 
formulation.  The  method  proposed  uses  the  Kinetic  Energy  Preserving  scheme  proposed  by  Jameson  and  a 
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Figure  15.  Lift  and  Drag  history  -  Sr  =  1.5 


modified  version  of  the  JST  artificial  dissipation  model  in  the  far  field  to  ensure  a  non  oscillatory  solution. 
The  resulting  code  proved  to  be  robust  and  extremely  low  dissipative.  Applications  with  coarser  grids 
containing  “only”  1024  x  256  cells  (results  not  presented  in  the  present  document)  showed  that  even  though 
the  far  field  results  were  largely  degraded,  the  time  history  curves  of  Cg  and  Cd  were  still  relatively  close 
to  the  ones  obtained  with  fine  grids.  On  an  other  hand  it  appears  that  the  code  can  be  easily  modified  to 
deal  with  airfoil  motions  more  complex,  like  a  combination  of  pitching  and  plunging  motion  obeying  non 
sinusoidal  variations.  These  remarks  lead  us  to  think  that  this  code  would  be  particularly  well  suited  in 
the  study  of  optimal  motions  of  rigid  airfoils  (in  a  sense  of  propulsion  efficiency)  at  relatively  low  reynolds 
numbers. 
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